TP GNS : ScaraBIS ou l’études des matrices de comptes

Les Librairies

La plupart des difficultés d’installation des librairies venaient de uwot, dont une version précédente était nécessaire et le couple SeuratWrappers et monocle3 dont leur installation demandent un token. Afin d’obtenir le token, je vous conseille de regarder ce guide https://gist.github.com/z3tt/3dab3535007acf108391649766409421 . Une fois le token bien enregistré dans .Renviron et R redémarré, la série suivante d’installations devrait permettre l’installation des packages problématiqes:

install.packages("devtools")
library(devtools)
devtools::install_version("uwot", version = "0.1.10", repos = "http://cran.us.r-project.org")
library(uwot)
install.packages("Seurat")
library(Seurat)

BiocManager::install("Rsamtools")
install.packages("R.utils")
library(Rsamtools)
library(R.utils)
remotes::install_github('satijalab/seurat-wrappers')
library(SeuratWrappers)
devtools::install_github('cole-trapnell-lab/leidenbase')
BiocManager::install("limma")
BiocManager::install("batchelor")
library(limma)
librairy(batchelor)
devtools::install_github('cole-trapnell-lab/monocle3')
library(monocle3)

Pour toutes les libraries restantes, un simple install.packages("...") suffira.

#devtools::install_version("uwot", version = "0.1.10", repos = "http://cran.us.r-project.org")
library(uwot)
library(Seurat)
library(tximport)
library(ggplot2)
library(corrplot)
library(network)
library(Signac)
library(SeuratWrappers)
library(monocle3)
library(Matrix)
library(patchwork)
library(gridExtra)
set.seed(1234)
setwd("~/mydatalocal/scarabi/results")

Récupération des données

Récupération des matrices des 7 échantillons:

sampsf <- c("SRR8257100","SRR8257101","SRR8257102","SRR8257103","SRR8257104","SRR8257105","SRR8257106")

files <- file.path(
  paste("~/mydatalocal/scarabi/data/quant/",sampsf,"/alevin/quants_mat.gz", sep=""))

txis <- lapply(files, function(f) tximport(files = f, type="alevin"))

Création des objects Seurat associés :

seu_objs <- lapply(seq_along(txis), function(i){
  # min.cells = 3, min.features = 200 from Ryu et al.
  s <- CreateSeuratObject(counts = txis[[i]]$counts , min.cells = 3, min.features = 200, project = sampsf[i]) 
  })

Pour l’instant, on le récupère que les samples du Wild Type (non mutants):

scarabWT <- merge(x = seu_objs[[1]], y = unlist(seu_objs[2:4], recursive = F), add.cell.ids = sampsf[1:4])

Etudes et pré-traitements de nos données

Ensuite, on étudie nos échantillons afin de ne conserver qu’une partie de nos données, en enlevant les données ne nous semblant pas significatives:

scarabWT[["percent.mt"]] <- PercentageFeatureSet(scarabWT, pattern = "ATM")
scarabWT[["percent.chloro"]] <- PercentageFeatureSet(scarabWT, pattern = "ATC")
quant <- quantile(scarabWT[["percent.mt"]]$percent.m,0.95 )
scarabWT <-subset(scarabWT, subset =  percent.mt < quant & percent.chloro < 0.2)

Voici le résultat obtenu :

VlnPlot(scarabWT, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.chloro"), ncol = 4)

plot1 <- FeatureScatter(scarabWT, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(scarabWT, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

Ensuite, on normalise nos données et l’on garde qu’un certain nombre nombre (pour ma version, 2000) des cellules les plus intéressantes:

scarabWT <- NormalizeData(scarabWT, normalization.method = "LogNormalize", scale.factor = 10000)
nb_cells = 2000
scarabWT <- FindVariableFeatures(scarabWT, selection.method = "vst", nfeatures = nb_cells)

Maintenant, on lance une ACP sur nos données :

all.genes <- rownames(scarabWT)
scarabWT <- ScaleData(scarabWT, features = all.genes)
scarabWT <- RunPCA(scarabWT, features = VariableFeatures(object = scarabWT))

Et l’on obtient ceci :

VizDimLoadings(scarabWT, dims = 1:2, reduction = "pca")

DimPlot(scarabWT, reduction = "pca")

DimHeatmap(scarabWT, dims = 1:15, cells = 500, balanced = TRUE)

L’étude des composantes principales a pour but de déterminer la dimension de notre problème. Afin de pouvoir déterminer le plus de types cellulaires possibles dans nos données, il est préférable d’avoir une dimension élevée. Cependant, chaque nouvelle composante principale doit aussi avoir un intérêt à discriminer nos données et doit donc être significative. Alors, pour se faire, on utilise deux analyses supplémentaires : JackStraw et l’ElbowPlot. Voici les graphiques obtenues :

scarabWT <- JackStraw(scarabWT, num.replicate = 100)
scarabWT <- ScoreJackStraw(scarabWT, dims = 1:20)

JackStrawPlot(scarabWT, dims = 1:15, xmax = 1, ymax = 1)

ElbowPlot(scarabWT)

Ce que l’on y obverse sur l’ElbowPlot est à partir de la 10ème composante principale, les suivantes ne sont plus si intéressantes que ça. On a donc décidé (arbitrairement) de se restreindre à une dimension de 10.

Génération des cartes UMAP

Afin de générer la carte UMAP (avec RunUMAP), on a besoin au préalable de convertir nos données en clusters, via un graphe pondéré en calculant les voisins de chaque cellules.

scarabWT <- FindNeighbors(scarabWT, dims = 1:10)
scarabWT <- FindClusters(scarabWT, resolution = 0.8)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11321
## Number of edges: 342018
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9284
## Number of communities: 21
## Elapsed time: 1 seconds

scarabWT <- RunUMAP(scarabWT, dims = 1:5, return.model = TRUE, umap.method = "uwot-learn")

On obtient alors les résultats suivants :

p8 <-DimPlot(scarabWT, reduction = "umap")
p8

ggsave(p8,file = "~/mydatalocal/scarabi/results/image/UMAP.png",width=30,height=15,units="cm")

Annotation de la carte UMAP

Afin d’annoter la carte UMAP en associant les différentes régions avec des types cellulaires, nous avons utilisé deux méthodes.

1ère méthode : via les marqueurs génétques

Le but de cette méthode est d’utiliser des marqueurs génétiques pour nos différents types cellulaires.

Markers = read.csv("../data/Markers.csv", sep = '\t', h=T)
Table = table(Markers$Preferential.expression.in.root)

Voici le tableau représant le nombre de marqueurs par type cellulaire :

Table

## 
##                                     Cortex 
##                                         12 
##                                 Endodermis 
##                                          9 
## Epidermis - Non-hair Cells (Atrichoblasts) 
##                                         11 
## Epidermis - Root Hair Cells (Trichoblasts) 
##                                         10 
##             Epidermis and Lateral Root Cap 
##                                          5 
##         Meristem (Young) Endodermis/Cortex 
##                                          7 
##                       Root Cap (Columella) 
##                                         17 
##              Stele (Young Vascular Tissue) 
##                                         13

Ensuite, pour chacun des types cellulaires, on affiche sur la carte UMAP seulement l’expression des gènes de ce type cellulaire.

Markers$Locus<-gsub(" ","",Markers$Locus)
Markers$Preferential.expression.in.root<-gsub("/"," ",Markers$Preferential.expression.in.root)

lm<-split(Markers,Markers$Preferential.expression.in.root)


system("mkdir -p ~/mydatalocal/scarabi/results/image")
system("rm -r ~/mydatalocal/scarabi/results/image/*")

output <- lapply(names(lm),function(x){f<-FeaturePlot(scarabWT, features = lm[[x]]$Locus)
  ggsave(f,file=paste0("~/mydatalocal/scarabi/results/image/",x,".png"),width=40,height=40,units="cm")
  })

datascore <- data.frame(lapply(names(lm),function(x){score=colMeans(scarabWT@assays$RNA[lm[[x]]$Locus,])
  }))
names(datascore)<-make.names(names(lm))
scarabWT <- AddMetaData(scarabWT, metadata = datascore)
g<-FeaturePlot(scarabWT, features=names(datascore))
g

ggsave(g,file = "~/mydatalocal/scarabi/results/image/type_cellulaire.png",width=40,height=40,units="cm")

2nd méthode : via des échantillons des régions associées

Via des prélévements on a aussi accès à l’expression des gènes par type cellulaire

samps<-read.table("~/mydatalocal/scarabi/data/Counts_Salmon/metadata_Li2016.txt", sep="\t")
#On trie les samples par types cellulaires
samps <- samps[order(samps$V3),]
#On supprime les types cellulaires non déterminants
samps <- samps[!samps$V3%in%c("whole root","whole root 1","cycloheximide treatment","cycloheximide mock"),]
ech<-samps$V1

Ensuite, pour chacun de ces samples on récupère les matrices de compte et via à une table, on réassocie ensemble les noms des gènes.

files <- file.path(
  paste("~/mydatalocal/scarabi/data/Counts_Salmon/",ech,"/quant.sf", sep=""))
files<-files[file.exists(files)]
tx2gene<-read.table("~/mydatalocal/scarabi/processed_data/txp2gene.tsv")
names(tx2gene)<-c("TXNAME","GENEID")
tx2gene<-unique(tx2gene)
head(tx2gene)

##          TXNAME    GENEID
## 1   AT1G01010.1 AT1G01010
## 7  AT1G01020_P2 AT1G01020
## 15 AT1G01020_P6 AT1G01020
## 21 AT1G01020_P1 AT1G01020
## 30 AT1G01020_P3 AT1G01020
## 38 AT1G01020_P4 AT1G01020

On peut alors récupérer les données et on obtient alors la matrice de compte (normalisée) suivante :

txis <- lapply(files, function(f) {
  tab<- tximport(files = f, type="salmon", tx2gene=tx2gene)
  return(tab$abundance)
  })
tabpur<-as.data.frame(txis)
ech2=sapply(files,function(f){strsplit(f,"/")[[1]][6]})
#Changer le nom des colonnes du tableau 
names(tabpur)<-make.names(ech2)
head(tabpur)

##           SRR3664382 SRR3664383 SRR3664385 SRR3664395 SRR3664406 SRR3664417
## AT1G01010   1.716811   1.896926   1.545711  15.843074  20.268393  27.714170
## AT1G01020   7.369509  11.111502   9.280125  17.836873  19.710903  14.096914
## AT1G01030   0.987394   1.474313   0.275222   3.689367   3.250175   3.115996
## AT1G01040   9.021629   8.665490   4.201520   5.220210   4.976830   6.818153
## AT1G01046   0.000000   0.000000   1.405640   0.000000   0.000000   0.000000
## AT1G01050  80.577198  96.215644  71.923622  80.756158  70.374419  67.141841
##           SRR3664396 SRR3664397 SRR3664398 SRR3664419 SRR3664420 SRR3664421
## AT1G01010   1.871991   3.618997   0.644235  11.358003  13.906525   8.793016
## AT1G01020  18.363256  16.313037   7.771199   9.363583   9.615792  10.181512
## AT1G01030   0.227622   0.047495   0.050757   0.094970   0.213155   0.238876
## AT1G01040   6.728216   5.175364   2.719236   5.278344   5.429668   5.547827
## AT1G01046   0.000000   0.000000   0.630508   0.000000   0.000000   0.819966
## AT1G01050  57.092115  71.099493  33.008388  92.716870  85.178594 103.502638
##           SRR3664402 SRR3664403 SRR3664404 SRR3664428 SRR3664435 SRR3664436
## AT1G01010  23.222836   5.580218   9.812438   3.259710   3.046452   4.292104
## AT1G01020   5.655989  10.831351   7.492238   3.077802   2.955736   1.872432
## AT1G01030   6.632473   3.191752   3.568764   0.050170   0.000000   0.234531
## AT1G01040   8.307703   7.156862   9.470113   2.188932   1.946090   1.110394
## AT1G01046   0.000000   0.000000   1.257230   0.000000   0.000000   0.000000
## AT1G01050  28.472910  46.433552  36.656115  39.044528  34.581970  41.300804
##           SRR3664405 SRR3664407 SRR3664408 SRR3664422 SRR3664423 SRR3664424
## AT1G01010   1.595829   1.669632   1.277080  30.262334  25.156696  24.216769
## AT1G01020   9.652126   6.374632   8.140137   7.261154   4.492351   6.229815
## AT1G01030   0.271967   0.308613   0.746134   0.150571   0.042560   0.035130
## AT1G01040   3.506512   3.004830   3.717348   5.564117   5.734301   5.665490
## AT1G01046   0.000000   0.548677   0.000000   0.000000   0.000000   0.000000
## AT1G01050  59.719122  50.988552  59.998486  48.434264  35.376794  36.978165
##           SRR3664374 SRR3664375 SRR3664437 SRR3664376 SRR3664377 SRR3664378
## AT1G01010  15.959187  18.168734  16.162023  57.228042  93.080597  83.976778
## AT1G01020   3.445357   3.176128   2.993439   8.075620  10.339301   7.353513
## AT1G01030   1.131589   1.439347   1.989967   3.381933   7.935872   5.471184
## AT1G01040   2.952577   2.444447   2.646486   9.455829  11.725951  11.487360
## AT1G01046   0.000000   0.000000   0.000000   1.121207   1.812935   6.681024
## AT1G01050  40.327258  25.524147  34.209081  30.594107  28.945098  25.546687
##           SRR3664389 SRR3664391 SRR3664425 SRR3664426 SRR3664427 SRR3664379
## AT1G01010   4.049538   9.052426   2.310489   1.541885   2.024963   4.263242
## AT1G01020   6.610719   3.517385  17.869868  16.397251  21.811314   5.658542
## AT1G01030   0.176140   0.000000   0.291288   0.322293   0.495996   0.358554
## AT1G01040  18.197588   7.864732   3.984504   3.071081   3.582128   2.855099
## AT1G01046  15.835439   0.614266   0.349925   0.000000   0.000000   0.000000
## AT1G01050  75.779156  59.369302 132.357691 101.977525  91.140975  57.979350
##           SRR3664380 SRR3664381 SRR3664372 SRR3664373 SRR3664384 SRR3664386
## AT1G01010   3.593588   5.582260  64.629978   9.211267  39.910340  55.078668
## AT1G01020   6.424684   5.219054   6.489138   2.905927  13.995652   6.056696
## AT1G01030   0.267133   0.488625   0.188868   0.658552   0.420462   0.000000
## AT1G01040   2.350481   2.728850  14.268800   6.399822  16.890305  15.307501
## AT1G01046   0.000000   0.000000   0.584967   0.000000   2.268211   0.000000
## AT1G01050  51.868707  48.501560  48.210757  35.352159  73.961593  43.024911
##           SRR3664387 SRR3664388 SRR3664392 SRR3664393 SRR3664394 SRR3664412
## AT1G01010  53.791810  37.618442   9.763332  18.649549   5.724260   0.604879
## AT1G01020   1.283277   5.794505  19.191938   8.081400   8.419938  17.998218
## AT1G01030   0.074136   0.000000   0.223327   0.069081   0.177317   0.622511
## AT1G01040  20.716409  17.283298   6.199051   5.724517   7.265130   4.821658
## AT1G01046   0.000000   0.721670   0.000000   0.748399   0.000000   0.000000
## AT1G01050  55.027919  43.526574  58.147394  42.592574  49.841970  48.721253
##           SRR3664413 SRR3664414
## AT1G01010   2.017441   2.845328
## AT1G01020  14.953029  16.290637
## AT1G01030   0.362785   1.603667
## AT1G01040   4.778282   8.002214
## AT1G01046   0.000000   0.000000
## AT1G01050  43.536514  50.942496

Enfin, pour chaque cluster, on regarde l’expression moyenne des gènes dans le cluster. Puis on ne garde que les gènes en commun, on calcule la matrice de corrélation entre les clusters et les samples et finalement on change le nom des samples (aux colonnes) par le type cellulaire associé. On obtient alors ceci :

avg.e <- AverageExpression(scarabWT)
scarabWT_avg=data.frame(avg.e)

genes_scarabi <- rownames(scarabWT_avg)
genes_li <- rownames(tabpur)
genes_common <- genes_scarabi[genes_scarabi%in%genes_li]

countsLi_norm_avg_alt_sum_c <- tabpur[genes_common,]
scarabWT_avg_c <- scarabWT_avg[genes_common,]

corLi_scarab_spearman <- cor(scarabWT_avg_c,countsLi_norm_avg_alt_sum_c,method="spearman")

colnames(corLi_scarab_spearman) <- lapply(colnames(corLi_scarab_spearman), function(name){samps[samps$V1==name,3]})

corrplot(corLi_scarab_spearman, method="color", is.corr=F, tl.col = as.color(colnames(corLi_scarab_spearman)) )

Alors, pour chaque cluster, en prennant par exemple le type cellulaire avec la plus grande expression, on peut annoter la carte UMAP :

cluster.id <- max.col(corLi_scarab_spearman)
cluster.id  <- colnames(corLi_scarab_spearman)[cluster.id]
names(cluster.id) <- levels(scarabWT)
scarabWT <- RenameIdents(scarabWT, cluster.id)
p7 <- DimPlot(scarabWT, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
p7

ggsave(p7,file = "~/mydatalocal/scarabi/results/image/UMAP_label.png",width=15,height=15,units="cm")

Etude du mutant

Afin de comparaitre nos mutants avec le Wild Type, on recommence toutes les opérations jusqu’à réobtenir la carte UMAP :

files <- file.path(  paste("~/mydatalocal/scarabi/data/quant/SRR8257105/alevin/quants_mat.gz", sep=""))
txi <- tximport(files = files, type="alevin")
  scarabmut <- CreateSeuratObject(counts = txi$counts , min.cells = 3, min.features = 200, project = sampsf[6]) 
scarabmut[["percent.mt"]] <- PercentageFeatureSet(scarabmut, pattern = "ATM")
scarabmut[["percent.chloro"]] <- PercentageFeatureSet(scarabmut, pattern = "ATC")
quant <- quantile(scarabmut[["percent.mt"]]$percent.m,0.95 )
scarabmut <-subset(scarabmut, subset =  percent.mt < quant & percent.chloro < 0.2)
scarabmut <- NormalizeData(scarabmut, normalization.method = "LogNormalize", scale.factor = 10000)
scarabmut <- FindVariableFeatures(scarabmut, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(scarabmut)
scarabmut <- ScaleData(scarabmut, features = all.genes)
scarabmut <- RunPCA(scarabmut, features = VariableFeatures(object = scarabmut))
scarabmut <- FindNeighbors(scarabmut, dims = 1:10)
scarabmut <- FindClusters(scarabmut, resolution = 0.8)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2371
## Number of edges: 67404
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8897
## Number of communities: 14
## Elapsed time: 0 seconds

scarabmut <- RunUMAP(scarabmut, dims = 1:5, return.model = TRUE, umap.method = "uwot-learn")
p6 <-DimPlot(scarabmut, reduction = "umap")
p6

ggsave(p6,file = "~/mydatalocal/scarabi/results/image/UMAP_mut.png",width=30,height=15,units="cm")
avg.e <- AverageExpression(scarabmut)
scarabmut_avg=data.frame(avg.e)
genes_scarabi <- rownames(scarabmut_avg)
genes_li <- rownames(tabpur)
genes_common <- genes_scarabi[genes_scarabi%in%genes_li]
countsLi_norm_avg_alt_sum_c <- tabpur[genes_common,]
scarabmut_avg_c <- scarabmut_avg[genes_common,]
corLi_scarab_spearman <- cor(scarabmut_avg_c,countsLi_norm_avg_alt_sum_c,method="spearman")
colnames(corLi_scarab_spearman) <- lapply(colnames(corLi_scarab_spearman), function(name){samps[samps$V1==name,3]})
corrplot(corLi_scarab_spearman, method="color", is.corr=F, tl.col = as.color(colnames(corLi_scarab_spearman)) )

cluster.id2 <- max.col(corLi_scarab_spearman)
cluster.id2  <- colnames(corLi_scarab_spearman)[cluster.id2]
names(cluster.id2) <- levels(scarabmut)
scarabmut <- RenameIdents(scarabmut, cluster.id2)
p5 <- DimPlot(scarabmut, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
p5

ggsave(p5,file = "~/mydatalocal/scarabi/results/image/UMAP_mut_Label.png",width=30,height=15,units="cm")

Comparaison entre le Wild Type et le mutant

Maintenant que l’on a les cartes UMAP des deux, un moyen intéressant de les comparer est de projeter les cellules et cluster du mutant sur la carte du Wild Type. On obtient alors ceci :

scarab.anchors <- FindTransferAnchors(reference = scarabWT, query = scarabmut,
                                      dims = 1:30, reference.reduction = "pca")

scarabmut <- MapQuery(anchorset = scarab.anchors, reference = scarabWT, query = scarabmut,
                       reference.reduction = "pca", reduction.model = "umap")

p1 <- DimPlot(scarabWT, reduction = "umap", group.by = "seurat_clusters", label = TRUE, label.size = 3,
              repel = TRUE) + NoLegend() + ggtitle("Reference annotations")
p2 <- DimPlot(scarabmut, reduction = "ref.umap", group.by = "seurat_clusters", label = TRUE,
              label.size = 3, repel = TRUE, split.by = "orig.ident") + NoLegend() + ggtitle("Query transferred labels")
p1 + p2

ggsave(p1+p2,file = "~/mydatalocal/scarabi/results/image/projection_mutant.png",width=30,height=15,units="cm")

Comme on peut l’observer, toute la partie poilue a disparu chez le mutant. Alors, une analyse supplémentaire est possible : on peut étudier le processus de différentiation des cellules de ce groupe via le PseudoTime.

atricho.cds <- scarabWT[, scarabWT$seurat_clusters %in% c(15,11,6,9,17)]
atricho.cds <- as.cell_data_set(atricho.cds)
atricho.cds <- cluster_cells(cds = atricho.cds, reduction_method = "UMAP")
atricho.cds <- learn_graph(atricho.cds, use_partition = TRUE)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%

# meristem to root
meristemcells <- names(scarabWT$seurat_clusters[scarabWT$seurat_clusters==17])
atricho.cds <- order_cells(atricho.cds, reduction_method = "UMAP", root_cells = meristemcells[22:22])
scarabWT <- AddMetaData(
  object = scarabWT,
  metadata = atricho.cds@principal_graph_aux@listData$UMAP$pseudotime,
  col.name = "Trichoblast"
)
p=FeaturePlot(scarabWT, c("Trichoblast"), pt.size = 0.1)
p3 <- plot_cells(
  cds = atricho.cds,
  color_cells_by = "pseudotime",
  label_branch_points=FALSE,
  show_trajectory_graph = TRUE,
)
p3

grid.arrange(p, p3, ncol=2)

ggsave(grid.arrange(p, p3, ncol=2),file="~/mydatalocal/scarabi/results/image/peudoTime.png",width = 10, height = 5)